{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "380f2a32-f419-43e1-91e5-b8e8fabac6b3",
   "metadata": {},
   "source": [
    "# 1. Bring in full data and precipitation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "2352af84-124c-4134-9dbd-7873d1c26ff0",
   "metadata": {},
   "outputs": [],
   "source": [
    "from basic_libraries import os, np, plt, xr\n",
    "from ape_funcs import get_full_data_fpaths, surf_temp_giver"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "2387bc04-cd85-46bb-a237-98e5dd84ef54",
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys\n",
    "sys.path.append('.')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "efac2c94-5a39-447b-9151-905802c9f307",
   "metadata": {},
   "outputs": [],
   "source": [
    "surf_temps = surf_temp_giver()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "74f3e2b6-1d56-461e-8d30-6f7742cec6f7",
   "metadata": {},
   "outputs": [],
   "source": [
    "#note these are 1D filepaths!\n",
    "root_dir = r'...'\n",
    "fname = '_hourly3d.nc'\n",
    "root_names = ['Ts280', 'Ts290', 'Ts300', 'Ts310']\n",
    "\n",
    "varis = ['prec_mp','w']  # variables: vertical velocity, cloud condensate, temperature\n",
    "fpaths = [os.path.join(root_dir, Ts + fname) for Ts in root_names]\n",
    "\n",
    "data_dic = {k: xr.open_dataset(fpath, decode_times=False)[varis] \n",
    "            for k,fpath in zip(surf_temps,fpaths)}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d8157722-83f3-43c5-aafc-89a612003897",
   "metadata": {},
   "source": [
    "### 1.0.1 Set up to use height, pressure and temperature as vertial coordinates "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "10d88d99-607d-4f1f-b043-eb6d0e324708",
   "metadata": {},
   "outputs": [],
   "source": [
    "oneDfpaths = os.path.join('...')\n",
    "\n",
    "pparte,sparte = 'Ts', '_warm_1d_mean.nc'\n",
    "fpaths = [os.path.join(oneDfpaths,pparte + Ts[:3] + sparte) for Ts in surf_temps]\n",
    "\n",
    "height_var,temp_var = 'zfull','tabs'\n",
    "\n",
    "oneD_dictionary = {k:xr.open_dataset(fpath,decode_times=False)[[height_var,temp_var]] \n",
    "                   for k,fpath in zip(surf_temps,fpaths)}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "7fc5d10a-a9c4-4b4b-870a-16cd88dd8aee",
   "metadata": {},
   "outputs": [],
   "source": [
    "#set up height and temperature as vertical coords\n",
    "temp_1d_dict = {k:oneD_dictionary[k][temp_var].mean(dim='time') for k in oneD_dictionary.keys()}\n",
    "zfull_1d_dict = {k:oneD_dictionary[k][height_var].mean(dim='time') for k in oneD_dictionary.keys()}\n",
    "\n",
    "from A_up_functions import extract_strictly_increasing_segment as esis\n",
    "#strictly increasing temperatures in the tropopause, first element is coordiantes, second is mask\n",
    "temp_vert_coords = {k:esis(temp_1d_dict[k].values) for k in temp_1d_dict.keys()}\n",
    "\n",
    "#assign temperature coordinates that are strictly increasing to data\n",
    "da_dic = {\n",
    "    k: data_dic[k].assign_coords(\n",
    "        zfull=(\"pfull\", zfull_1d_dict[k].values),\n",
    "        temp_profile=(\"pfull\", temp_1d_dict[k].values)\n",
    "    )\n",
    "    for k in data_dic\n",
    "}\n",
    "\n",
    "\n",
    "#mask Data to be in strictly increasing temperature portion of tropopause\n",
    "trop_ds_dic= {\n",
    "    k: da_dic[k].isel(pfull=slice(int(np.argmax(temp_vert_coords[k][1])), None))\n",
    "    for k in da_dic\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ec6f2a8c-1a58-407a-a8dc-fd412707d09a",
   "metadata": {},
   "source": [
    "# 2. Make mask using data from aloft: w and maybe condensates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "a55c58dc-6111-4043-bcf3-3d5382e45ded",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Processing datasets: 100%|████████████████████████| 4/4 [00:02<00:00,  1.55it/s]\n"
     ]
    }
   ],
   "source": [
    "from tqdm import tqdm\n",
    "condensates = False\n",
    "\n",
    "Ts_names = ['280K', '290K', '300K', '310K']\n",
    "\n",
    "\n",
    "w_umbrales = [(-0.15, 0.7) for i in range(4)]\n",
    "\n",
    "# Threshold dictionaries\n",
    "w_thresholds = {k: val for k, val in zip(trop_ds_dic.keys(), w_umbrales)}\n",
    "\n",
    "mask_holder_up = {}\n",
    "mask_holder_down = {}\n",
    "\n",
    "for Ts_name, (k, dataset) in tqdm(zip(Ts_names, trop_ds_dic.items()), total=len(Ts_names), desc=\"Processing datasets\"):\n",
    "\n",
    "    w = dataset['w'].values\n",
    "    condition_w_up = w >= w_thresholds[k][1]\n",
    "    condition_w_down = w <= w_thresholds[k][0]\n",
    "\n",
    "\n",
    "    joint_up = condition_w_up\n",
    "    joint_down = condition_w_down\n",
    "\n",
    "    mask_up = np.where(joint_up, 1, 0)\n",
    "    mask_down = np.where(joint_down, 1, 0)\n",
    "\n",
    "    da_up = xr.DataArray(mask_up, coords=dataset['w'].coords, dims=dataset['w'].dims)\n",
    "    da_down = xr.DataArray(mask_down, coords=dataset['w'].coords, dims=dataset['w'].dims)\n",
    "\n",
    "    mask_holder_up[k] = da_up\n",
    "    mask_holder_down[k] = da_down"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "44f5915d-7150-42f6-aff9-357549008215",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Filling Holes: 100%|█████████████████████████████| 4/4 [00:00<00:00, 21.78key/s]\n",
      "Filling Holes: 100%|█████████████████████████████| 4/4 [00:00<00:00, 21.33key/s]\n"
     ]
    }
   ],
   "source": [
    "# Collapse to 2D\n",
    "mask_2d_up = {k: mask_holder_up[k].any(dim='pfull').astype(int) for k in mask_holder_up.keys()}\n",
    "mask_2d_down = {k: mask_holder_down[k].any(dim='pfull').astype(int) for k in mask_holder_down.keys()}\n",
    "\n",
    "mask_2d_upF = fill_holes(mask_2d_up)\n",
    "mask_2d_downF = fill_holes(mask_2d_down)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "47ba13b1-9324-4734-9242-1fb217aa78cd",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Filling Holes: 100%|█████████████████████████████| 4/4 [00:00<00:00, 20.85key/s]\n"
     ]
    }
   ],
   "source": [
    "mask_2d_union = {\n",
    "    k: (mask_2d_upF[k] + mask_2d_downF[k]).clip(max=1)\n",
    "    for k in mask_2d_up.keys()\n",
    "}\n",
    "\n",
    "mask_2d_unionF = fill_holes(mask_2d_union)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "dda8a817-9cfb-4808-ae71-db62b3dabb2a",
   "metadata": {},
   "outputs": [],
   "source": [
    "categorical_mask_2d = {}\n",
    "\n",
    "for k in mask_2d_upF.keys():\n",
    "    up = mask_2d_upF[k]\n",
    "    down = mask_2d_downF[k]\n",
    "    \n",
    "    # Create new categorical mask\n",
    "    combined = up + 2 * down  # up=1, down=2, both=3, neither=0\n",
    "    categorical_mask_2d[k] = combined"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "486f39a1-4861-43c5-b62c-3e41ee57dec7",
   "metadata": {},
   "outputs": [],
   "source": [
    "lpf_mask = {k: xr.where(categorical_mask_2d[k] == 0, np.nan, 1) for k in categorical_mask_2d}\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "9c3d72ae-df49-4e38-b7ce-75291d2ea7f9",
   "metadata": {},
   "outputs": [],
   "source": [
    "#revaporation data, will be used to mask precip based on detection of mass flux event aloft\n",
    "\n",
    "prec_array_dic = {k: da_dic[k]['prec_mp'].values for k in surf_temps} #these are the precipitation fields for hourly averages, and re-evap on"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "823c3f2e-7068-4c90-8202-d01192a3e095",
   "metadata": {},
   "outputs": [],
   "source": [
    "#this is where the mask is made\n",
    "lpf_masked_prec = {k:prec_array_dic[k]*lpf_mask[k] for k in lpf_mask}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "754ca454-3ddd-49dd-8983-4398157946e4",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "10.509\n",
      "23.192\n",
      "39.744\n",
      "63.447\n"
     ]
    }
   ],
   "source": [
    "for k in lpf_masked_prec:\n",
    "    print(np.round(np.nanmean(lpf_masked_prec[k].values),3))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "556d8424-06b2-4503-ac77-16f237043c2b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.82\n",
      "2.69\n",
      "3.67\n",
      "4.7\n"
     ]
    }
   ],
   "source": [
    "for k in prec_array_dic:\n",
    "    print(np.round(np.mean(prec_array_dic[k]), 2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "777ca353-b870-47bb-b301-59e9b55ba2a9",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "280K: NaN=0.868, Non-NaN=0.132, Sum=1.0\n",
      "290K: NaN=0.908, Non-NaN=0.092, Sum=1.0\n",
      "300K: NaN=0.926, Non-NaN=0.074, Sum=1.0\n",
      "310K: NaN=0.939, Non-NaN=0.061, Sum=1.0\n"
     ]
    }
   ],
   "source": [
    "for temp in ['280K', '290K', '300K', '310K']:\n",
    "    nan_frac = np.isnan(lpf_mask[temp].values).mean()\n",
    "    print(f\"{temp}: NaN={nan_frac:.3f}, Non-NaN={1-nan_frac:.3f}, Sum={nan_frac + (1-nan_frac)}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1ca85bdd-b70f-4dff-bca0-b0c33b72b98a",
   "metadata": {},
   "source": [
    "## Make plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "d0d250f9-3647-4467-b2f3-c3da599cc7c4",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib\n",
    "import matplotlib.cm as cm\n",
    "import matplotlib.colors as mcolors\n",
    "from matplotlib.colors import LogNorm\n",
    "\n",
    "import cmocean"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "25c5ec2b-0000-4e34-89de-6ac0ad96cc91",
   "metadata": {},
   "outputs": [],
   "source": [
    "def make_masked_precip_panels(fig=None, positions=None, precip_data_dict=None, \n",
    "                             temp_cases=None, time_steps=None):\n",
    "    \"\"\"\n",
    "    Create two masked precipitation panels with shared colorbar.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    fig : matplotlib.figure.Figure, optional\n",
    "        Figure to plot into. If None, creates new figure for preview.\n",
    "    positions : list, optional  \n",
    "        List of two GridSpec positions for the panels. If None, creates subplots.\n",
    "    precip_data_dict : dict\n",
    "        Dictionary of DataArrays keyed by surface temperature.\n",
    "    temp_cases : list of str\n",
    "        Two temperature cases to plot, e.g. ['280K', '290K'].\n",
    "    time_steps : list of int\n",
    "        Two time step indices, one for each temperature case.\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    dict\n",
    "        Contains 'axes' (list of two axes) and 'cbar_ax' (colorbar axis).\n",
    "    \"\"\"\n",
    "    # Preview mode\n",
    "    if fig is None:\n",
    "        fig, axes = plt.subplots(1, 2, figsize=(8, 3))\n",
    "        ax1, ax2 = axes\n",
    "        preview_mode = True\n",
    "    else:\n",
    "        ax1 = fig.add_subplot(positions[0])\n",
    "        ax2 = fig.add_subplot(positions[1])\n",
    "        preview_mode = False\n",
    "    \n",
    "    # Extract and plot data\n",
    "    norm = LogNorm(vmin=1, vmax=2000)\n",
    "    \n",
    "    for i, (temp_case, time_step) in enumerate(zip(temp_cases, time_steps)):\n",
    "        ax = ax1 if i == 0 else ax2\n",
    "        \n",
    "        # Get data for this temperature case and time step\n",
    "        da = precip_data_dict[temp_case].isel(time=time_step)\n",
    "        \n",
    "        # Plot with masked regions\n",
    "        im = ax.imshow(da.values, cmap=cmocean.cm.rain, norm=norm,\n",
    "                      origin='lower', interpolation='nearest')\n",
    "        \n",
    "        # Set title with temperature\n",
    "        ax.set_title(f\"$T_s = {temp_case[:-1]}\\,K$\", fontsize=9)\n",
    "        ax.tick_params(bottom=False, left=False, labelbottom=False, labelleft=False)\n",
    "    \n",
    "    # Add shared colorbar to the right of the second panel\n",
    "    if not preview_mode:\n",
    "        # Create colorbar axes to the right of ax2\n",
    "        divider = make_axes_locatable(ax2)\n",
    "        cax = divider.append_axes(\"right\", size=\"5%\", pad=0.01)\n",
    "    else:\n",
    "        # For preview, create colorbar on the figure\n",
    "        cax = fig.add_axes([0.88, 0.11, 0.02, 0.77])\n",
    "    \n",
    "    cbar = fig.colorbar(im, cax=cax)\n",
    "    cbar.set_label('mm/day', fontsize=9)\n",
    "    cbar.set_ticks([1, 10, 100, 1000,2000])\n",
    "    cbar.set_ticklabels(['1', '10', '100', '1000','2000'])\n",
    "    \n",
    "    if preview_mode:\n",
    "        plt.show()\n",
    "    \n",
    "    return {'axes': [ax1, ax2], 'cbar_ax': cax}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "a30ddcd3-d26e-415b-8730-661632713d01",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAqUAAAEWCAYAAABfbu50AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAABObUlEQVR4nO3dd3wUdf4/8NfupvcESIMAQWoUATEUQaUjAoLKIYqKoiDcAoeoWKginifqqZwodvT7O/RsFFFBTk4RiVRRQERKpAgJJb23/f0x2c/skk0yM1tmd/N6Ph77yOQzM5/57GYz+9n3pxksFosFREREREQ6MupdACIiIiIiVkqJiIiISHeslBIRERGR7lgpJSIiIiLdsVJKRERERLpjpZSIiIiIdMdKKRERERHpjpVSIiIiItIdK6VEREREpDtWSomIiIhId6yU6uSDDz5AREQEIiIiEBYWBoPBIH6PiIjAxx9/7LJrlZeXY8qUKUhNTUVkZCQ6d+6Md955x+6YP//8E2PHjkWzZs3QvHlzjB8/HufPnxf7KysrMWPGDMTGxiIuLg4zZ85EVVVVnWtVVlYiJCQEW7ZsEWmff/454uPj8corr7jsORERefI+CgAzZ85ESkoKoqKi0LJlS8yePRsVFRVi/yuvvIKrr74awcHBGDt2bJ3zeR8laoSFdPfBBx9YWrZs6bb8i4qKLAsWLLAcPXrUUlNTY8nIyLDExMRYNm3aJI4ZM2aMZcyYMZbCwkJLQUGBZfTo0ZYJEyaI/QsXLrR069bNcubMGcuZM2cs3bp1szz55JN1rvXTTz9ZDAaDJS8vz1JTU2NZtGiRJTk52bJt2za3PT8iInffRy0Wi+XXX3+1FBUVWSwWi+X8+fOWAQMGWJ566imx/9NPP7WsWbPGYjabLWPGjKlzPu+jRA1jpNQL/PTTT7jqqqvcln94eDiWLFmCyy67DAaDAX369MHAgQOxbds2cczx48cxfvx4REREIDIyErfddhv2798v9r/zzjuYP38+kpKSkJSUhHnz5uHtt9+uc629e/eiffv2sFgsGD16NL7++mvs3r0b/fr1c9vzIyJy930UALp06YLw8HAAgMVigdFoxJEjR8T+W265BWPHjkXz5s0dns/7KFHDWCn1Anv37kWPHj0UH//Xv/4VMTEx9T5sK5uOlJWVYefOnbjyyitF2pw5c/Dxxx8jPz8feXl5+OCDDzB69GgAQG5uLk6fPo3u3buL47t3746TJ08iPz+/znOJiIhAeno6UlJS8O233yIpKUnxcyMi0sJT99F//OMfiIiIQHx8PH7++WfMnDlT0fV4HyVqHCulXkDtN/xXX30VeXl59T769+9f77kWiwX3338/OnTogFtuuUWk9+vXD+fOnRN9nXJzc/H4448DAIqKigAAMTEx4njrdmFhoV3+e/fuxYkTJ2AymbB8+XIEBQUpfl5ERFp56j762GOPoaioCL/++iumTZuGxMRERdfjfZSocayU6uzkyZO4ePGi25udAKlC+te//hWHDx/G2rVrYTRKf/6amhoMHToU/fr1Q1FREYqKitCvXz8MGzYMABAREQEAdt/mrduRkZEirbq6Gj///DNWrVqFiIgIPPDAA25/TkREnryPWnXp0gXdunXDPffco+h43keJGsdKqc727t2L5s2bIyUlRaRVV1fjgQcewIABA9CvX786zUjTpk2zG2F66eP777+vcx2LxQKz2YwdO3bg66+/RnR0tNiXk5ODEydOYNasWQgLC0NYWBhmzpyJHTt24MKFC4iNjUWrVq2wb98+cc6+ffuQkpJil89vv/2GkpIS9O/fH+vWrcPGjRvx7LPPuvDVIiKqy1P30UtVVlba9SltCO+jRI0L0LsATZ2jflC//PILCgsL8e233wKQIpm2Vq5ciZUrV6q6zowZM/DDDz9gy5YtiI2NtdvXvHlztG/fHitWrMCiRYsAACtWrECrVq1Eh/17770XTz/9tOho//e//x33339/nefStm1bxMbGIjY2FuvWrcOAAQPQqVMnh9OjEBG5gifuo0VFRfj4449x8803Izo6GgcOHMDSpUsxfPhwcUxVVZV41NTUoKysDEajUTS/8z5K1DBGSnXmqB9U165d0b59e/zlL3/B66+/LprZtTpx4gReffVVHD58GG3atBGRgGnTpolj1q1bh71796Jly5ZISkrCzp07sX79erF/wYIF6Nu3L7p06YIuXbqgX79+eOKJJ+yus3fvXrvnkp6ejrfffht33nkn9u7d69RzICKqjyfuowaDAatXr8Zll12GyMhIjBkzBiNHjsRLL70kjlm6dClCQ0Px9NNP4/PPP0doaKjoBgXwPkrUGIPFYrHoXQiyV1JSgrCwMFRXV6NTp044evSo3kUiIvIpvI8S+R4233uhyZMn49SpUygvLxcj4ImISDneR4l8DyOlRERERKQ79iklIiIiIt2xUkpEREREumOllIiIiMiLPPPMM0hPT0dkZCTi4+MxduxYHD582O6YsrIymM1mNGvWDBEREbj11luRnZ1td8zJkycxcuRIhIWFIT4+Ho888giqqqrsjvn2229x1VVXITg4GO3bt8eqVavc/fTqxUopERERkRf57rvvYDab8eOPP2Lz5s2orKzEsGHDUFxcLI558MEH8fnnn+Pjjz/Gd999hzNnztgtH15dXY2RI0eioqIC27dvx3vvvYdVq1Zh4cKF4pjMzEyMHDkSAwcOxL59+zB79mzcf//92LRpk0efrxUHOhERERF5QHl5OcrLy+3SgoODERwc3OB558+fR3x8PL777jtcd911yM/PR4sWLbB69WqMGzcOgLQaWJcuXZCRkYE+ffrgq6++wqhRo3DmzBkkJCQAkBaNePTRR3H+/HkEBQXh0UcfxRdffIEDBw6Ia02YMAF5eXnYuHGji5994xRNCVVTU4MzZ84gMjISBoPB3WUiIoUsFgsKCwuRnJzs9OTg5F68jxJ5JzX30bKyMlRUVGi+1rJly/D000/bpS1atAiLFy9u8Lz8/HwAQFxcHABgz549qKysxJAhQ8QxnTt3RuvWrUWlNCMjA127dhUVUgAYPnw4pk+fjoMHD6JHjx7IyMiwy8N6zOzZszU/R2coqpSeOXPGbk1hIvIup06dQqtWrfQuBjWA91Ei79bYfbSsrAwR4VGorqnUfI3ExERkZ2cjJCREpDUWJa2pqcHs2bPRr18/XHHFFQCArKwsBAUFISYmxu7YhIQEZGVliWNsK6TW/dZ9DR1TUFCA0tJShIaGqn+STlBUKY2MjAQA9Fp4L45XBIn0Y48vc0+pvMDD78lLv225KHUKvqNluLTvtkUeK8eeLW8DABLaXCfSWl3WwWPX90W5OYUAgNi4SKfzsHImL3cqKChASkqK+B8l72X9G506dQpRUVE6l8Z/vfifJQCAB29b2MiRvufcmSyxffzgWgBAn6HT6jnanvWeNvb/PSXSDh+qBgBYaoODbdqbxL4PB98MAIhtdYVI89b7oLOU3kcrKipQXVOJ1kk9YTSaGjzWkZqaapw8uwchISGq7gFmsxkHDhzAtm3bVF/T1yiqlFqbmjY/9JLdC/nNl/8EAFzVZ4pI8/U3rfUfd3Z3uRKYcHwPAGBKv8kAgOoquenNHc931lsPi+344EAAwO2GUpHGD7SGueL18bXXmM3B3s/6N4qKivK595cvWTTleb2L4Da275v2neeqOtf6ufXh4FEi7avLdgAADhRIg2f6N4sR+9qkXQMAyM8+bnP9lo1eZ+OG58R2i+BoAEDPoVNVlVUvSu+jRpMJRqOGBTE13KZnzJiBDRs2YOvWrXZR3MTERFRUVCAvL88uWpqdnY3ExERxzM6dO+3ys47Otz3m0hH72dnZiIqK8niUFODoeyIiIiLlDAbtD4UsFgtmzJiBNWvWYMuWLUhNTbXb37NnTwQGBuKbb74RaYcPH8bJkyfRt29fAEDfvn2xf/9+nDt3ThyzefNmREVFIS0tTRxjm4f1GGsenqahqk9ERETUNBkMBk2tU2rOMZvNWL16NdatW4fIyEjRBzQ6OhqhoaGIjo7Gfffdhzlz5iAuLg5RUVGYOXMm+vbtiz59+gAAhg0bhrS0NNx1111YtmwZsrKyMH/+fJjNZtGPddq0aXjllVcwd+5cTJ48GVu2bMFHH32EL774QvXzcwWnKqWDb5wDADh1VJ7QNTauEwDg0XelpoXmQXIf1EcmLnXmch4hmuMvHyHSZnUYCABokdjCI2X4PFse2dcsWNqef1d/l+Vv21/SXd0tXNGvU0n+jujZhcQTry0RkVriftTpepE2yWb7Uqd/WQcA2J9zQqSFRDYHABRlS5/5FkuN2FeQKzXzHysuEmk1l0w46ei+bXufzDwodZVLvbxn/U/EG6iMetqdp9Brr70GABgwYIBd+rvvvot77rkHAPDiiy/CaDTi1ltvRXl5OYYPH45XX31VHGsymbBhwwZMnz4dffv2RXh4OCZNmoQlS5aIY1JTU/HFF1/gwQcfxMsvv4xWrVrhrbfewvDhw9U/PxdgpJSIiIhIKQM09Q9Vc46SKeRDQkKwYsUKrFixot5j2rRpgy+//LLBfAYMGICffvpJeeHcyCWV0pT2ncT25NcfBABU1b6e4Sb1I9TcQW0Ey/4Yz0a8roqR37kVNXX3a41CNhRdJNdgdJTIP/2+Q5pI/PdzvwAAogPlKX2uvWGWLmXSwvYe1dBngnUgz4UKeaL39dvfAgAcKJDOyyqXl6s01n5sBRnlz6+ssjIAQOq2zwAAiWlDxb6XN0hzdb57Uh7Ee31zaZjLS4nyDDNxzbxwUKAHIqVNFSOlREREREqxUuo2rJQSERERKWQwGGAwahjoZGGltDEuqZRmHtgttt954EUAwEsfLAAAtLBZqeD7jcsB6NPU4ahZtbFO155mLc+n5pfr3Qc432zviefo7mtY82eXBCLyhI69bwAADHzqKwDAnwvq3qd93Z/7Pxfbm/78FQCw5YJ8jy2olvrlHciVKlclZfK5wbVjmlMibfucSQN1M3OPAQDCzsnN8v86JjXbl52UuwCszpKa71uGyIOi5960wK6MXtFFipFSt2GklIiIiEgxD4x0aqJcUilNveLqOmmzb3/KwZGeZ42kvfblMyLtQrm0bu30rlKn6w69hnm+YA409A3QmW+HaqOK174oRbI7hMmD1KwRcG/T0OvibZFwIvJdfV+Q7outw6TfF773mNi3ZNI/PFKGf66Wo4Zz7lD3GXvs5+0AgMu6XSPSLmbuAgDsPJEBAPgy+7zY91WWFBUtyJYjn4baEeGG2oHMpip5hHhJhPR5URAiV7xOl0nLmK7POgsAyK1cJ/bd31ZqRV0uj6MCyqX8MkvkaRFP7vsUALDx7BEAwKMTnxb7/u9TKaJ6163z6zzfOW8/AgD4533P1dnnNNZJ3YaRUiIiIiKl2HzvNn5ZKbWNkFm/CW65KE/o2ypY+kZXVV5Q59wZbz0EAHjl/hfcWURdKI0Sfv/gcjeXxHlKpsVyRVTU3YsAEJH3sl0Ypn+cNAWUdW3uB0fP83h5lEZHHbUSWSOkR3Z+LdJ+PiNNVr89JwcAsCdfjorm1348tjgvT9kUWCFFPiuDpM9Qg81Umjkm6fW5UCS3sJVXSZWwI0XSeT/kystdtq6NqN7VTj7eOp1U92h5Gqhfck8DAP5zRnpOT8//m9hXsLT+fr15VVX17nOaxhWdWCltnF9WSomIiIjcgpFSt2GllIiIiEgp9il1G6cqpd7atGlfnnQAwH97DhIp69c/CwBIvLzu2q5BBqlxZvmHC0XarAlL6hxH+rp08Ja73oPe9t4mIs+xXa3wufbLAACDXpIGPHnzvaGhslVXyfM4FVRJg34rahept53MKSBIqkFZbObjNFZJzfAhtYOFTeXyKKVqk/TZeTFEngYyp9iamZTHBZtBUEUR0jXLLXIz+xUR0rxSRTZN783CpPyGNreuoGUzD1UD3Do410OR0q1bt+K5557Dnj17cPbsWaxZswZjx44V+y0WCxYtWoQ333wTeXl56NevH1577TV06CBPvZWTk4OZM2fi888/h9FoxK233oqXX34ZERER4phffvkFZrMZu3btQosWLTBz5kzMnTtX/fNzAWPjhxARERERAMBo1P5Qobi4GN26dat3bftly5Zh+fLlWLlyJXbs2IHw8HAMHz4cZWVyxX3ixIk4ePAgNm/ejA0bNmDr1q2YOnWq2F9QUIBhw4ahTZs22LNnD5577jksXrwYb7zxhrbXxkkGi8ViaeyggoICREdHY8pyMy4a5U7Jt7VMAAAMvW6mSLv0G5rSSd/dPX0PpwfyLHdEMOetelTOv1L6Jv3qFP0GpDmzoIErHN2zBYVFxbhqwE3Iz89HVJQXrhFNgvU+yr8V6cn2vpV1cBMAYMuZ/QCANVm5Yt/Z2nrN4dNyFaFZlhQZDSsoAQCYaiOmAJCfEAMAqAiSK17G2ghsaZjUKFsVLEcKw2Kk49pHy/lfFSUdlxYZLtL6NpPqGS3irwQAlBWcFvviOw6RfibF1/+EFVD6v2k9rl3XITCZ1Dc0V1dX4fj+/2q6BxgMBrtIqcViQXJyMh566CE8/PDDAID8/HwkJCRg1apVmDBhAg4dOoS0tDTs2rULV18tTd25ceNG3HjjjTh9+jSSk5Px2muvYd68ecjKykJQkBSpfuyxx7B27Vr89ttvqp+jsxgpJSIiIlLKALkJX9VDOr2goMDuUV5e3uDlHMnMzERWVhaGDBki0qKjo9G7d29kZEjzzmZkZCAmJkZUSAFgyJAhMBqN2LFjhzjmuuuuExVSABg+fDgOHz6M3Fz5S4qnsFJKREREpJSmCqncDzUlJQXR0dHi8cwzzzRywbqysrIAAAkJCXbpCQkJYl9WVhbi4+2jyAEBAYiLi7M7xlEettfwJFXx5y0Xq3CmXG6+L66WCtzt8HciLbbvKLtz3N2s+fuPX4rtFh2vrfeabKp3jpLmeHc3Z1fY9DRpHhTo8vydoba7giu6N3x86BuUlar/hk1ETZfdPad2sO9go1QV+KNku9hlrB1QlB0nx67yy6RoWkixdN8xGuQBScFlUlN+sM0tqTxYuk+HlEoDpCpq5LxKAqUKWqk83gaxQVI5rm3RUqS17zmxbrn15uRAp1OnTtk13wcHB9d3RpPDSCkRERGRUk5GSqOiouweWiqliYmJAIDs7Gy79OzsbLEvMTER586ds9tfVVWFnJwcu2Mc5WF7DU9SFSm1WICKYjladbBAeoH3Z+8XafE51wNw/K3m1LHfAQApl3Wss0/rt6A3D30rtp/rc6OmPPydp6bucnf+swdOF9st26a69VpKNPZ8//wjEwAwYY00NckfxfI366c6NwMAjBn0oOL8AOChdx4R2/mV1ahgpJSInBQc3QoAEB0gt4TGBUoxq9YR8kRRZ1Ok/WcCpShfVK480Mk6wKnGZgqp6tqpoIzVUr2h2mQzvVTtpSJtGr26REph08R28hSOVl41WNnJeUrT09NhMplgNpthNps1FSE1NRWJiYn45ptv0L17dwBSX9UdO3Zg+nTps7Jv377Iy8vDnj170LNnTwDAli1bUFNTg969e4tj5s2bh8rKSgQGSn+MzZs3o1OnToiNjdVUNmdw8nwiIiIipZxsvt+1a5ei0fdFRUU4evSo+D0zMxP79u1DXFwcWrdujdmzZ2Pp0qXo0KEDUlNTsWDBAiQnJ4sR+l26dMENN9yAKVOmYOXKlaisrMSMGTMwYcIEJCcnAwDuuOMOPPnkk7jvvvvw6KOP4sCBA3j55Zfx4otunOe1AaoqpVdGGhAUKv9+6DFp3dmft6wSaQ19cwmPaTwU7Gjd+vY2E99v+uJ5AMCunIsAgD358nxcqz97GgBwxy2eX5PYFdwV0XRFft7QnycsqrneRVDFGs29KV66+RwsLBb72oVHa8ozzCT3uGkZEoIyDdOSEBEB8n29KEf62S5c7uCZUyH1F82plFtjalDbUpoo3YfyouVmZ2sQ1HbJ+eoK+xknTUFyRa5F7UdK/xh51Hf3mCQAQERsUr1l9Q4aK6Uqw6u7d+/GwIEDxe9z5swBAEyaNAmrVq3C3LlzUVxcjKlTpyIvLw/9+/fHxo0bERISIs7597//jRkzZmDw4MFi8vzly5eL/dHR0fj6669hNpvRs2dPNG/eHAsXLrSby9ST+IlGREREpJTBKD20nAflzfcDBgxAQ1PJGwwGLFmyBEuW1L/qZFxcHFavXt1gsa688kp8//33jRTeM1gpJSIiIlLKQ833TZGqSumIhGa4IVg+5bcf1gIAYlumKzo/rln9fwRHnZitzfaj//U3kZYWITUX3NyyjVSGoiNiX06Ffwz60HulIG/lq6/FIxOXAgBOHj4o0lp3ulxTXk9NelZsH9r2GYqKS+CbnVWIyFuktO8EAIg5JDelR9QOekoOlgc/GQ3S1E6VtSs1GR1UzGpsAnsBtYOe2obXHg95Z8sQKWrYKy5GpKVeNaFOfl5533dyoBPVj5FSIiIiIqUMBliciJRS/VRVSq9KSENEeJj4vXO/sS4riPXbkKOI6eczXxbbyz9cCADoM0zqh1H4hbz2+dCRD7msPESuFtmitUvz69L/FhQUFLg0TyJqeqyfu22iU0RaVnkpACDcZpqolNpRTIlBFdIxFdViX7CxboUrIUg6t18zaWqhSbcuqHPMhXMXxbZXRkUdcbL53hVTQvkrRkqJiIiIlGKfUrdhpZSIiIhIKScrpVQ/VZXSjn1Gqa7dqx20Y3uMo3k77xr2iN3x/tRk7zNNF05y9Hf11KpTevLn50ZEvst6b4q9Th5oZB2KufUruftcbqU0L/hTkx4FALywer7YF2iUBi61DJUnM08JjQEAdLj6zjrXtN7zTQFBdfZ5PVZK3YaRUiIiIiKFLAbpoeU8apjbK6XORMMcHcdok/s5GmxmpfX1t83zxL5PAABvnJWXT6uunSB4wpU3AQDaXdlH03WUlsMX3kecGoyI9Na192Sxfel96KE7ljZ4rpL7rU/e24xG6aHlPHCgU0MYKSUiIiJSyKJxSigLBzo1ipVSIiIiIqXYp9RtPFopdWWYXs9mWEfN297cBPHr9x8BANKuHa85D2efX1HuWbG9P/dPAMC2nCKR9muh1HzfLGgjAOABNzXfW/lC07i3louIfM+uzSsBAOlDp6k6z5n7kN8129din1L3YaSUiIiISClGSt3GayqlWiOfvhDx0pvaCKk7ItpZx7eItFCTtMpHWkSwSKu0SFONXN+yu8uu7QjfI0TUFBVUluhdBP9hgMZKqctL4nc0DB8jIiIiapqsA520PABp9H1aWhpWrFih8zPxPl4TKVUbwfKWiJe3lMNbWV+fwspSkTZu7OMAgB77tom0yuILAIDO/caqyn/Lly8CAAbd+KAzxSRqkr796iUAwIARs3UtByljvd8B6u95g2+c4+riNFmW2oeW8wCOvm+I11RKiYiIiLyfxj6lbL9vFCulRERERAo5O08p1c9rKqXns84DAFokthBpmQd2AwBSr7hac77umL6JTfbqOWpquqx7f7fk60urNjkjN6cQBQX1r75F1Bg22/uWMF9cJ14nr3+0GAAQESBXcybeMt8leVuMBliMGiqlGs5parymUkpERETk7RgpdR9VldLzWefd1jnXNkJqZY2QHrMZEOOK6Br5j6YSFXXk3OHvUFTMaV5I9tHavwMAWodFi7Q+w5rm2tpbv3pZbF834m86lkTm7P2qqf4t1XjjoycBAMdr7409YqIbOlwbA7R1D609Jz09HSaTCWazGWYz/6a2GCklIiIiUsjZSClH39ePlVIiIiIihdh87z6qKqWOmtg9Ia51N7HdlJtrqS5H74Om8t4wBYXBWKFltjzyV2Em6ZZ+rqxI55Loz5NN9o4G1Co5vqncqzxp6vhFAIAP1zwNAOjf7TbXX8TJ5nuqHyOlRERERAoxUuo+PlEpdebbJL+Jkr9q33MQCgoK9C4GuYgrBnSOGj3XVcUhJxXnS9McFmf9KtJqqisAAIlpQwE0PmUhI6raTbh5ntvythikh5bzqGE+USklIiIi8gaMlLoPK6Uusm/Lu2I7Or4rACAmuRMA13zLtU71AgB7cvMAAM/euwyA/bdtrdc6nXlcbLdKbdfo8e5YlMAZ329cLrarLTUAvGdicG97rcg7cbo737fhf/I0VGtrF4QJdTBh+uV/7AEATLvxsQbzO7L7/wEA/nf+NADg0YlPu6Sc5BxOnu8+rJQSERERKcTme/cx6l0AIiIiIl9hbb7X8gCkyfPT0tKwYsUKnZ+J9/FopNQfmzHf+lhaPeJIUbFIKz1+EAAwb8RDtSnOP8fxY58Q24//fRYA4PdXpClP1sx42eE5aihpsrflbX+3+Gi5/OVlOTqWpH7e9poRkXa2n2eFF6Xm9X8cPSfSDv9Zu2GSp20LCZUqJd+ESFN2xWx+XuybMPThOvlGhDYDADQLkvMl/TkbKeXk+fVj8z0RERGRUhorpZyntHEuqZQqHWjjrZGioz9tFds15dI32I59blR07okSaX3dL8+XibSzxdI7r+P3rwEAZkxY4nQZ/+/TpWLb2lf6XO3E6Qe3fij2XX7dBKev5Ys69R2ldxHq5a3veyJS7tKWvpqaarH9TsZ7AIDf/5CjolH5lQCAyPxSkVYZHAgAONc6DACw7Giu2Ncx8k0AwFV9poi0pMtHAAAGx6pryfJH3jQ9lsVgADj63i0YKSUiIiJSyAJoinpy/b3GsVJKREREpBAjpe7jkkqpN4TTnVFw8bDYvmrIlAaOlBza9pnYvq55PADgbNmfIu1imNSs0y2muauKiFED5XWcTQZpTs6oQKkpKPmKkYryOHdW6ixfXiw3GVUUSXPp6T1H4qVNY77+niIi/1ZZJg9u/TFP6sYVUFkj0qIvSF3BAnPzRVpQaDAAoEVt3eRMULjYt/z4CQDA6OJ/irThV90OAEi94mpXFt0nOfpMuPp5adDv7oeX19nnThYDtEVKWSdtFCOlRERERArVGAwwMFLqFk5VShvqeOxNnZIb0+ryMaqO79L/FrGdWPs8Xfk91tHAMdvXccSAWXbHK32N45OkqO6PX38s0i6US9/2i3PlFZ2uHHi3pjISkTL8//F9a354Q2zvvihVNqoCbab+rq1/WGqjowBQFiP9rfNipbQuLeRehrEBJumYankAVdbxLQCAiNgk6Ri+V+xMbd1Ml+syUuo+jJQSERERKcRKqfuoqpRu3vgySkzyqxpd26fx2v5/rXOsL32js0YQtXDH82wsT2ev2WeYWWxv+fLFOvv1iHI7ey1fWJiB0TGy4t/fewx6SWp5Crf5bPt8Zt0FSax/M+v/cXZ5udhnnabPEizncbZ1DAAgtKRKpBU2CwIAdEiWfn82rbXY175lHwDAiTO75LSeE+2uTfamjl+ky3U50Ml9GCklIiIiUoiRUvcxNn4IEREREQHyMqNaHgCQnp6OtLQ0rFixQt8n4oVURUqH3vA3fLBRnq7ifG3zRcs9H4i0PblnAQDrsqT1x1sEyfXeF8dLKxuxKcJ7DLrxQU3nNfY39FQXAOtqXBnH5VW54kNCAAC9+j4g0vR8z/V9QWoevC4uRKQ9e+8yvYpDRDasTe+OmuwbMvX66WK7RfA7AICvz+eJtF8Lpc++zhHyQKe7UhIAAEOvm1lvvnGtu6kqB+lAY6TUes6uXbsQFRXlyhL5DTbfExERESnEPqXuo7pS2ixI/tb3wPjHAQB7NstTY+RWVgAAzpZLU10YDfJkwrl/Hqrd6iLSGopg+dK0UmTP3X+zwxkbAACZF6X31MoT58S+UGlmFTwX8pFIO2WUBuUpme7KVf65egEAILz2v6y0uqaBo4lID9PbJqk63tG97Y6hDwMAbjwjL8Ty+4lvAQDtErqLtObteqvKl7wT+5S6DyOlRERERApZjNA0IsfCUTyNYqWUiIiISCkn+5RS/VRXSseNfbxOWruet4vtO2oHOrX5RVo1KNQkXyI6qaOqa7E5wzGt3Rr8qTvEsQtSs/3/LpwHABwtlP/b44KlriP5FfK8oN163+fB0knm3PGU9LP2d0dzqRKRvnq0G6zpvOK8bLEdGtUcABCT3Emk9bLZJj+jsU+ppnOaGAaTiYiIiEh3Lmm+t428WbeHx85q8LimwDrQxRoxcxWtr6MnX/9Lo4KuvvaNox8BAGS89ygA4BqbJZBTQqRVUzqm3ezSazqrqb3/iXxBu669NJ3Xql17sX3yiNRyE9mslaJzeS/wbRzo5D7sU0pERESkFPuUuo3bKqX8Juj6CKk3OLpnCwAgtFmqSGvZNrXOcZ76+/eNk/pyPTXpEY9cj4joUq07dKmT5qgPeVP9XHzjoycB6LdWvVuwgukWjJQSERERKcVIqduwUkpERESkkNbB96yUNk51pTT7z7NiO6GlupUwvIWzUyO5u1nGNn9va+6prigBAJRcOCYnOmi+9xTrgCcianqO7P6v2O5w9RDdyuHoM8W6bXs/96dp+dTwq2Z7oDZSqmVKKJeXxO8wUkpERESkkMEoPVTjJJyNUl0p9dXoqC1XfkttKt94rd/w4ztdD6DpPW9bTeW5E3k7PaOjthxFRRs6jnwc+5S6DSOlRERERApprZNS41gpJSIiIlKKkVK38boeDrk5heKh1c6vXxOPjE2vIGPTKy4sodQEY334EmdfVyIiqp/tZ4OjB/kH6+h7LQ81tm7ditGjRyM5ORkGgwFr166122+xWLBw4UIkJSUhNDQUQ4YMwZEjR+yOycnJwcSJExEVFYWYmBjcd999KCoqcvIVcB+vq5QSEREReStPVUqLi4vRrVs3rFixwuH+ZcuWYfny5Vi5ciV27NiB8PBwDB8+HGVlZeKYiRMn4uDBg9i8eTM2bNiArVu3YurUqc48fbfymuZ769Qe7xz4WqQ9c88yTXn1GjZdbFvXJCb1nexto6oXju8AAEQkdJYSmsi3fkev2enjRwHYr32th+P7d6LQi7/xEpF3W/L+4wCAhXc/o3NJGrf5ixcAAIcK80XarAlLdCmLp+YpHTFiBEaMGOFwn8ViwUsvvYT58+djzJgxAID3338fCQkJWLt2LSZMmIBDhw5h48aN2LVrF66++moAwL/+9S/ceOONeP7555GcnKzhSbgXI6VERERECjkbKS0oKLB7lJeXqy5DZmYmsrKyMGSIPANFdHQ0evfujYyMDABARkYGYmJiRIUUAIYMGQKj0YgdO3Y49yK4ie6RUms07ufTuwAA23LlsPPUN+YAAN6Y+k/N+Ttak9iX6DnZsv01eztVDn+aWik8JkG3a586elhs//nnjyguLmvgaCKi+vlChNRq6MiHpJ86lwNwPlKakpJil7xo0SIsXrxYVVZZWVkAgIQE+8+jhIQEsS8rKwvx8fF2+wMCAhAXFyeO8Ta6V0qJiIiIfIVJ4+T5ltpzTp06haioKJEeHBzsopL5PjbfExERESlkcOIBAFFRUXYPLZXSxMREAEB2drZdenZ2ttiXmJiIc+fO2e2vqqpCTk6OOMYVXBl19ZpI6bixUofrcQ72ObMWvDesNezNa9kr5avldgc9X4uU9p3EdlV5EQc6ERF5mFFj872l9pz09HSYTCaYzWaYzWZNZUhNTUViYiK++eYbdO/eHYDUV3XHjh2YPl0a7N23b1/k5eVhz5496NmzJwBgy5YtqKmpQe/evTVd15G2bdti6NChmDx5MkaPHo2AAO1VS6+plBIRERF5O2f7lO7atcuu+b4+RUVFOHr0qPg9MzMT+/btQ1xcHFq3bo3Zs2dj6dKl6NChA1JTU7FgwQIkJydj7NixAIAuXbrghhtuwJQpU7By5UpUVlZixowZmDBhgktH3h87dgz/93//hyeeeAIPPPAAJk6ciMmTJ6Nr166q89KlUuqKyGdD59seY7FYVJbOu7giKqfkNXM3vSOtqz97GgBwxy3znM7LG6LvAJB6eU8UFBToWgYif/HDxn+J7cjgaABAeJw87dtl3a7xeJnIOxkMUrRUrRqV5+zevRsDBw4Uv8+ZIw3+njRpElatWoW5c+eiuLgYU6dORV5eHvr374+NGzciJCREnPPvf/8bM2bMwODBg2E0GnHrrbdi+fLl6gvfgJYtW+Kxxx7DY489hoyMDLz33nu4/vrr0a5dO0yePBl33XUXIiOVfV6yTykRERGRQgZonBKq9vz09HSkpaXVOym+1YABA2CxWOo8Vq1aJZXDYMCSJUuQlZWFsrIy/Pe//0XHjh3t8oiLi8Pq1atRWFiI/Px8vPPOO4iIiHD9i1IrOTkZSUlJiI6ORl5eHj799FOkpKTg/fffV3Q+m++JiIiIFDJqjJSqbb73FSUlJfjkk0/w7rvvYvfu3bj55pvx7rvvYsCAAQCk+VJHjhyJu+++u9G8PFopPfbzdgBAUETzOvsaagq13ad27XaDpo4frqW2mdfVzcLW/JryuveuaLa3cvR6av2bWVdUAXxrzkAib7Ppi+cBAMNHPqzp/H43zBTbR3ZKKwuyyZ4cMYLNzLYSExPRuXNnTJ48GWvXrkV0dLTd/r59+yoeWMVIKREREZFCzkZK/c22bdtw5ZVXNnjMV199pSgvl1dKD2dsAADUVEkrzSRePlzsc8W3TiURKb0HoLjbpQNt3vr4SbHv/r8s0qVMTZEr3meny+Tl5ca88jcAwHPp0lrHHXvf4HT+RE3F2dISAK4Z1Nih1zCXlIn8k8kIGDWESq0T7rtiSihvYlshLSwstBtgrrabAiOlRERERAqxT6m9EydOYMqUKdi2bRvKy8vt9lVXV6vKy+WV0k59R7k6S6qHNWLaM66VouP9PYLckCO7/wsAqCzNAQCkXTu+zjF6LHIwqXWqfP3KUgBAiw79PHJtIn9yz7iFeheBmgg239szm82IjY3F9u3bcf3112Pr1q1YvHgxRo4cqTovRkqJiIiIFOJAJ3sZGRk4ceIEIiIiYDAY0K1bN7z11lu49tprcf/996vKi5VSIiIiIoW0Tp5v8dNIqclkQnBwMACpD2lOTg6io6Nx6tQp1Xn5RKX09Y8Wi+0Hxi+u97imqm2Puk3RZN8cbwyQVrhIunxEvcfr0b3BdhoaIiLyfiaD9FCt9hx/G+jUvXt3bNmyBcOHDxfR0bCwMKSlpanOyycqpURERETeQGufUoufDnR6++23UVNTAwB4+eWX8fjjj6OgoADvvfee6rx8olJaqnD01sYNzwEAiqqqRNq4sY/Xd7jf8NUBTO4eWJRzar/YjmvdzW3XISKipkNrn1JL44f4pJSUFLHdvHlzvPnmm5rz8olKKREREZE3cDZS6g+UrmWvZGlRW6yUEhERESkUaABMBvVxT03TSHmpF1980e73gwcPIiwsDElJSTh79ixKSkpwxRVX+GeldPbtTyk67lhxEQAgwGa9+0tXP/IH3raWvTe8xravRebeDwEASZ3kQU3+9PcnIiL9MFIK/PTTT2J7/vz5GDlyJBYvXoyAgABUVlZiyZIlmvLlVFtEREREClmnhFL7MNiMvk9LS8OKFSv0fSIu8vrrr4sKKQAEBgZi4cKFWLlypeq8fCJSqpT5tifrpLkymugNEUFb3hKZ1MqV5bfNq7BNfwBAUoqyla48RY8Vo4iIyLW0DnSynuNvo+9DQ0Oxd+9e9OrVS6Tt27cPISEhqvPyq0opERERkTtJ85Sqb4v3p+Z7W3PnzsWwYcNw5513ok2bNjhx4gT+/e9/46mnlHW9tMXmeyIiIiKFtDTda+2H6gtmzJiBzz77DGVlZfjf//6HsrIyfPLJJ5gxY4bqvFRFSnNzCn0u5OyOJmI2w3q3yGZSs72jLgbu+nvNfvthAEBZ7QTCALBgqLRaU8u2qW69NhF5h/c+lSND+ZWVAIBZE7QN+CDv5Wzzvb9Yt24dhg0bhtDQUAwaNAiDBg1yOk9/e42IiIiI3IaRUsmKFSuQnJyM0aNH480338TZs2edzlNVpNQdkR5vGzzkyMnDBwEAH+xeDQBoFhQk9t3/l0W6lElv3vz38tSUWbb5J4UEAwCKq+TVx2qqKlx2rTMnTgAAQiPjRJo3/w2o6Tl55JDYPnLkKwDA4Bvn6FUcXUy6dYHeRSAPkPqUqj/P2qc0PT0dJpMJZrMZZrPZtYXzoK+//hqFhYXYtGkT1q9fj3nz5qFt27a46aabMHr0aHTr1k11nhzoRERERKQQR9/LIiMjMW7cOIwbNw41NTXYvn071q9fjwkTJqC0tBSjR4/GtGnTcPnllyvKT/dKqS9Ee44e+xoA8P5pKTLWXA6UYuhv0vrqbTp39Xi5CNj8xQsAgKEjH6qzz93vLdv87xkgdeguyf1TpKW07+Sya1kjpL7w/0JNU+sOXRxuE/kbkwEI0FAr9dfR91ZGoxH9+/dH//79sWzZMhw9ehTr16/HoUOHfKdSSkREROQrTLUPtWoaP8RnVVdX48SJEygqKrJLnzNHXRceVkqJiIiIFDIaDDBqmKdUyzm+YO3atZgyZQouXrxol24wGFBdXV3PWY6xUqrAoBsfBAD0ODEbABAfFCj2RcW3BcBpohrz7Vcvie0BI2bX2X/pgDelr6ejZntPsS1jQsskacP608X4niJfErPgbwCAvKde1rkkRK7HKaHszZo1C88++yxuv/12hIaGOpUXK6VERERECmmd3snfpoSyKisrwz333AOj0flqt+pKqT9GBDMP7BbbqVdcXe9xd6a0BACcKSsVaQd2vgsAuKLXvW4qnX/4o7igbtqhn8V2WHQCAODDNcsBAFcl9ZAP7HgtAO9+v/nC1GZE7nSidtAnAMzrzP8D8l+MlNozm8149dVXNa3gdClGSomIiIgUYp9Se1OmTMH111+PZ555BgkJCXb79u7dqyovVkqJiIiIFHK2+d5fJs+3GjduHFJTU3HzzTcjLCzMqbxUV0p9tXny9x+/FNvLD3wDADhXLo0KG9wiRuwbnyzNLenoeR4syAcAfJYlN0VHBlwAADxmeVukORrI407vfCKvrZwW1QwA0GeYd73R24VH10kLCA4X2x9//xoA4JXMXOn4U5vEvvfa9Qbg2bXslfDV/wXyP5u+eB4A0DxImpC759CpHi+D7VzNj3DeZvJjBmhrirfWY/1p8nwA+OWXX5Cbm4vAwMDGD26Ev3ZxICIiInK5AINB88Mf9enTB8eOHXNJXk2m+f738wfE9nc5lQCAIukH0iLLGjz35y2rpJ8FxQCAw/nyG6t1hEXKy4XrnKtVWSNPyXsgX4rc9tGrMPW4bsTf6qSFx8h9T2IDpWWyImu/aFXDIvYd+ekDAEDbtLEiLaQ28soBRkTA8JEP610E8gML3ntUbD816VkdS+LdTAbpoeU8f9SzZ08MHz4cd955Z50+pbNmzVKVV5OplBIRERE5y2AwwKAh6qnlHF+wc+dOtGvXDtu3b7dLNxgMrJTWZ9TouWL7w9OzAQBnyqUIY4+YWLHPUcQtq1SKPo5JigcAvD/tcTmvNU8DAPr1m+7aAqvwwPjFul3bGbav9TVdRgEAzpV/AgAorZajv5+dOQkAKD75okjrGCF1ph7RVupvGtt3lHsLS0QEYOfXUv/3XsP0u+fbGvSS/KG/ZfZyp/JidFQZTgll73//+x8AoKqqCiUlJU7l5a+vEREREZHLWUffa3n4ox9//BHdu3dHSEgIYmNj7R5qNZlIKREREZGzjDDACA3zlGo4xxdMmjQJt99+Oz744APPTwnlbkv/T24a35UvhYHXzXDt+sl/qV2f/Hy5NMDpKpsBNI5Y30i3jnm8zr4JN89zadmaqrZdugEAJiW0AwBkH/pG7PvyxB4AwA8lNqtCFUnvjRsscjM/EbnOkvfl+93Cu5/RsSTeJdAUpHcR7DjbZE/qcZlRe+fOncOiRYtc0meWzfdEREREChlqI6VqHwY/jZTecccdWL9+vUvy8miktKHJz6375t/l/m/k6d0mAACS27RRdPzQkQ+5sziKZGx6RWz3HS6tL+ttk8m7UmVlkdgODzABALpGhoq01HBp4v3O19zk2YIR+bm+L0gDZzIeYgTOkR6D79O7CHZOZx4X261S2+lYkqaDkVJ7S5cuRZ8+fbBs2bI6U0J99tlnqvLyuuZ7IiIiIm+ldSL8aj+dEurOO+9EUFAQ+vfv7399SomIiIi8lbOR0vT0dJhMJpjNZpjN3rUkuBbfffcdzp49i8hI51tqPVopbahpuaF9ts3UrmieVtps7w22fiUN8nrjj0yRtvZdac7VkYktAQDJCd1szugJwHeb8a3lzmveWaTF558GAETbrKvbJbqlZwtG1ESkR3vXQB5qGJvsPc/Z0fe7du1CVFSUq4ulm7S0NBQWFvpepZSIiIjIlxkM0kPLef7olltuwahRozB9+vQ6fUpvuknduA+fqJTqHfXTc331FhFSRDDA8IdIyyypAAC8clyKng4ryhP7rq+S9qF9X5GmttzesJ586hVXO9zWi6uj9d7wGhM5svz+5wG4/j3v7b7fKA/sOlsmDbQcP/YJvYpDXozzlNp7/fXXAQB///vf7dINBoN/VkqJiIiIvAFH39vLzMxs/CCFdK+UWr+NH9/zAQAgv7JY7OvYZQwA/fvM6Bkl6HLtOADAvNAYkfaf2onl3zstvVY/5OSIfTfnbwIATCu9INJir79D1TWbQlREKUfTbmmNcjrKi2jD58sAAKNGz9W5JJKm+j49VSI/72ATp/Cm+hk0Rkr9dZ5SV9K9UkpERETkKxgpdR9WSomIiIgUMhgMmpbUdMUynP7O5ZXSS5t+lDZxvnb8NwDA2fJqkTa9+hMAQKtUfZu1ju/fCQBo17WXbmVo3q632B5VcAoAsDb7JwDAhVL5uKMllR4tl79zZVcGdosgR7yl2d6qqb5P77hlnuZz9215FwDQfdC9rioOabRm3T8AANUWCwBg3NjHXX6NAGibPD+AzfeNYqSUiIiISCGOvncfl1dKL13L3hHbfVm/bgYA7C2oAgBcLJf/aCdKiuscr+RbvKunMtEzQmpl9zy61w5+KswGAPxeJD/fG5I6AACCgmM8VjZqWEP/C001KkXkT5pqhPT3HRvFdsfeN+hYEtnNYx4DAHz71UtuuwbnKXUfRkqJiIiIFGKk1H1YKSUiIiJSiKPv3cdtlVJrs+Sv338k0opKLwIAeg2bLh+YNhQAMO3MfgBARY080OmOoQ/b5dWYpjK/nvX1GHCtWfrpYJ9ap48fFdvhMQlO5eUKR3b/V2xb52VtFiSvyf2XgbMAAM3jm3m2YCpY34/nj/wAAMjJPyb2JSRbB651EmlsyiciX+ItTfaODBgx2215Gw0GGDW0xWs5p6lhpJSIiIhIIUPtQ8t51DC3V0rTrh2v6Lip4xe57JreHHFy5aAXVz7P0gvHxXZZ/pna/K9zWf5qff/HdrH9SZY04C0+SF7tKzHjTQByp/Yl78vTfrQOCwUA3DNuodvLeSnbv+/pX9YBABYe2AUAKKy2iH035UiD1EZWy5GG2Lhr6s33cMYGAECnvqNcV1giIlLNoDFSynlKG8dIKREREZFCxtqHlvOoYbpXSpva5ORFOVIUMjAkQqQltmqpV3GEDr2G6V0EOz3jWovtK87nAABiAk0irXlwmN3xC+9+xjMFU2F/zgkAwM5ca4r8LbljWG3U95A8pUqPs7sBAFf0kqaX+XDz82Lf1+fzAABrGCklUqTDM7PE9pHHl+tYEvI3XNHJfVhxJyIiIlIowAAEGAwaHuqus3XrVowePRrJyckwGAxYu3at3X6LxYKFCxciKSkJoaGhGDJkCI4cOeK6J6oDVkqJiIiIFLLOU6rloUZxcTG6deuGFStWONy/bNkyLF++HCtXrsSOHTsQHh6O4cOHo6yszBVPUxe6N983NSntOzV+EKHboHvE9v8b5N5rORp8prUriO15IwZIzYc/5z8FAMgurxT78qqkqc+2XSwQaaXVUlr43g8BAIFG+TtjlTxGiogUaB7U+DGkvxX/kQc5m297UseSKOepFZ1GjBiBESNGONxnsVjw0ksvYf78+RgzZgwA4P3330dCQgLWrl2LCRMmqC+gF2CklIiIiEgh6zylWh4AUFBQYPcoLy9XXYbMzExkZWVhyJAhIi06Ohq9e/dGRkaGy56rp3lNpNSV0SpXy9j0CgAgJjwZANCl/y2Kzju2bxsA4M+ze0TadSP+5uLSNW7zFy8AkCfuHXzjHI+Xwd1ybCKOcc2iVJ3rrveZNd+Z108FABw8tFbs25QlTQlVYrNYRKhJGsiVWyH9Lwxod63YlxK21y1lJPJXr/XoIbaz/zwLAEhomaRXcageCcEhehdBNWdH36ekpNilL1q0CIsXL1aVV1ZWFgAgISHBLj0hIUHs80VeUyklIiIi8nbOjr4/deoUoqLk4ElwcLDLyubr2HxPREREpJBB4yAnQ+1Ap6ioKLuHlkppYmIiACA7O9suPTs7W+zzRV4TKfWWpnqrTV/Ic0SeKJbmlOxc29TaRWEel3XvDwA4+uePqq5t25XBFa9LZKDUPFJY6bsj8hrjzfO/tWrXvvbnwyKt8vNlAIAzpaUiLdgkfUesrH2fBYTGin29+j7g9nIS+ZPug+7VuwikwLixj9dJyzwod3lLvbynJ4ujiNEgPbScBwDp6ekwmUwwm80wm82aypCamorExER888036N69OwCpr+qOHTswffp0TXl6A6+plBIRERF5Oy3TO1nPA4Bdu3bZNd/Xp6ioCEePHhW/Z2ZmYt++fYiLi0Pr1q0xe/ZsLF26FB06dEBqaioWLFiA5ORkjB07VnXZvAUrpfWorKkR29YBKHFhLQCoj2QOH/lwo8fYWrb+KbF9WXg4AODWwfLgJItFmh+oolQqR2MrQvUZpu2bmCP7trwLADhfmiPSUpPTAQDte1znsuuo5W2R9sZ0SroKABB74ZBIK6+pAgAkRkp/z7ZpPeqeSETk586c2i62vTFSajIaYNIQKlV7zu7duzFw4EDx+5w5Uj1g0qRJWLVqFebOnYvi4mJMnToVeXl56N+/PzZu3IiQEN8bPGbFSikRERGRQs5GSpU23w8YMEAEoRwxGAxYsmQJlixZoros3kpVpfT4/h3o3m+ou8riVfr1q9snw1PRuGfuWdbgfmukNjhM3dRH1vOceR77c/8EANRA/kfpHOW7naoB4NzZcwCA4osnRVrqFVe79ZodrpbmlusAeY65Myel6wcE1D/rt/UYACg6+yuKavs7ExH5i343zBTbG2r737dv1hkA0Pmam3Qpky1nJ89X2nzfFDFSSkRERKSQs5FSqh8rpUREREQK2a7OpPY8apiqSmm7rr3dVQ6v4y0DZxw1ubtiXXatfiuUymPbfN/x6NcAgJTLOjqdvyu6GKgVnxQvbVh/KuTqVchCI2Lt8rBdEzrYKE0X1T4iWqSZDEYUF/vvNF9ERKNGz9W7CHU4u6KTK6aE8leMlBIREREp5OyKTuxTWj9VldLn//MkxrXvJn6/cuDdLi9QU7T3v28CAA4VSCszTLxlvtjnLRFbq4gA6S1TbrNme5Wlpr7DNXP1AgLuZi2jM+W+9HjzbU82ePystx5GRWm5qmsQEWn16bpnAAC3jqk74X1TYtDYp9TAPqWNYqSUiIiISCH2KXUfVkqJiIiIFDIZDDBpqGBqOaepUVUpffi2Rfh5+yrx+8kj0mo0rTsoXQ2+aWqsSddgkLo/51VWeqxMWt1yWV8AQElJtkhr22M8AN9rcneW7XN86J1HAADtwsNEmrX53V2Dt7rHRKE0iAOdiMix1Z89LbYDalcTGj/2Cc35Dbp2htNl8gfORko50Kl+jJQSERERKeTsPKUc6FQ/VZXS3JxCtO1wg/jdFVMAuUP2n2fFdkLLJB1LImksQtZj8H3ST08Uxkmd+o6qk2aNBG76brlI259fAAAYlyK9R6zPsTG+GmEtrZYGe50vd93Ao8YizzcPehAFBQWYgWdcdk0i8h933DLPpfn56v3Z1QzQOPqeA50axUgpERERkUJc0cl9VFVKz/32DdJv8P5poIJCI1Qd7+q+kOfOZAEAAkPCXZant7F9zUoKLgAAPvjznEg7UCBNrl9W8xsA10eB9Zhkv74yAMCrU16od7/SMlqPL86T+useP/yF2Neu00gAwJ4Dn4q0iIBAFJewTykRkScZDNA4T6kbCuNnGCklIiIiUshoMMJoUL+mk5Zzmhq+QkREREQKWZvvtTwAafR9WloaVqxYofMz8T6qIqWdrhlrN8WEkk7UekwT5Mx1XNEsbG2292eOXp9rYuVuEwlBpQCA9NgYAMC5s3LTfrzKNea9le1r4Oh9o+Q9ZPv/UVNdBQC4f4M0YGz3Rbmt596sNwAA2eXytGHvT3sRBQUFAFw7mIGIiBqgcZlRcJnRRrH5noiIiEghk8EEk8Gk6TxqmOpKqW10dPmHCwEAUQFyNh0iYgEA/W6YCcD3Bvm4ory+9pydZX2+j058us4+dw1I8rbXWGt58s7+Lra/2L8OALA/X/o2XWwzu9TxEumXSJPc4+bulQ+iotR1U1AREbnDrLceBgAsv/95nUviGgaDUSx6o/Y8ahgjpUREREQKcUoo92GllIiIiEghRkrdR3WldMn7j4vt8hppFZsDFUUirbS6GgAQtHklACB96DSnCugJ3tYU7I20zrnpqddWjwF1rrBy53/E9p4Cac7RuGDp9+7R8rfqPjHSc6qBRaQ9OvFpFBQU4D9zXvNASYmItNGz2f7nLavEdrdB97gkT2crpenp6TCZTDCbzTCbzS4pk79gpJSIiIhIIWcrpRx9Xz/VldJWYaFie+UJaSWfI/lyRKekOh8AEBUYCABId6p45C3URh89Ha30peiore158opMp0ukn63CpJ+dI4LEvnYR0jRjMYHBIm3Nun+ghCs6EZEXcUWr1SPvzhXb5dVSi6zWaKuroqO2DDDAoGGadwP7lDaKkVIiIiIihdin1H1YKSUiIiJSiJVS91FdKU1v3k5sf3RGar63WWQGJ8qkUHtWmdSs+O1XL4l9A0bMrpOfpwfEkPvYNtvkZx8HAASESH/X8OgWYh//1rL4ILk550K5NIipTYh043ph8nNin/X/6NL/IWlFpyfdWkYi65zUBZXyzX7+Xc/oVRzSme293sp6X1d6f3/jI+m+9WuhnFdOpTRQum+s3N9y+m1eeH8zGKWHlvOoQYyUEhERESlkNAbAaFRffdJyTlOjfqDTlWPE9qLyPADA60HH5Axr13ZtE1Y7MCM42pnykY+KTmjX+EFNmDXS8HT3gSJtT/YBAIDJwZrKjloZiDxl1oQleheBvMDvP34JAJi1cyMA4JzNgnLvpfcCAHQdcKeivH7MzQMAvH+oWk6sllqLTl+WK5Kmay1sA+54bbbYXj39JfUZMFLqNqy2ExERESnEPqXuo7pSattfpG3aWADADSXviTTrhPrRgdJ0Nt0H3Sv2TXtzDgCgV2yMSJs8bqHaIpCbPffv+QCARyYuVXWe2r6iDfVL8maHtn0GAAiJbSPSUi/vqSoP8Ty7DBZpE/uNdbpsRESuZHufvmObFCE9vaMUABBQIUc5J2EnAOAbm9bUhu7n7cOl6SVNQcUiLTBAaiXqGG5yttgN0hQdteVkpJST59ePkVIiIiIihTh5vvuwUkpERESkkAEGbZVSTp7fKKcqpUkprQAA6TnyYI3Lul1T7/HHS6oAAJEBRc5cltysoWZ7R03uVt6+6pMzbJ/3a79tBwBcqNwq0ubnSIP9ki4fIdKUPD9Hx1w4dxEA0Dy+mbbCEvmh1z9aLLYfGL+43uPIfYy1daqKQKl5vSg8UOy7rVmwo1Pq9cSdfwcAXL7+WZEWFyQ16V97wyxniul+HOjkNoyUEhERESmlsfmeldLGqa6Ufr9xudi2fptxFB21Rpask6gDwF2t4gEA7SNi1F6WnKB2gQIlEyOvs/l226PLTbX7ujhVTm9m+9ptulABAMgulptiLgvbAwCYYxMp1YoRUqK6aiwWvYvQ5H1ykzRB0/vx7wMAukbHiH3X9f8rAPUtYGNuetQ1hfMkRkrdhpFSIiIiIoUMRhMMRvUzBGg5p6lhpZSIiIhIKUZK3UZ1pfRUidy0a10P2dFqH9YQfnFeuEiLDwkBAMSEJaq9LDnBlQOQhrwsddkIsBlEOMPwOQAgslkrzdfUk7W7QknBBZHWsm2qw2MA4PBjUheWeavkZqd+zZq7s4hETZ5XroHeBNjfy6V7vHnkE/oUxktw8nz3YaSUiIiISClGSt1GdaW0nc0gpZyc7EaPf+P7t8T2zvwSAMDGWQ+rvSx5iZ5RIXXSruw4EoBvRkcB4PBuqdN+RXWVSLtwPBoAEBgYBgAIa9ZB7Cu8KL0Gt7ZqL9KCAqUWAV96DYiI1Lj0/rZm3T/E9s1jHvN0cXRjMBpgMGqIlBo5T2ljGCklIiIiUoqRUrdRXSntM0xep7WPguPPVVSI7ZxK6ecrH8rr3d94+SgAQLuuvQAAf/z6k9hXlnsCANDssn4irUViC7VF9gjbyJv12+Tof/0NAPD5zJfdci09onLP3rvM49d0hy3fvyK2V53MAgBkV8hTzoxNiAAAdIqUXuPgvD/Evl/y8wAAf5SUibROEVJE9QoXlE3Pvy8RkVLuio72ekGePD8lRIoufmp27eeoM5ztU5qeng6TyQSz2Qyz2dzIWU0LI6VERERESjkZKd21axeioqJcXCj/wEopERERkUIcfe8+bq+Urpzyzzpph7Z9JratzfZWbdN6iO2DWw8D8N4me6DhplZXN9tbeWuzrqMuDN5mz+Y3AACbz50XaV/9XgMAMFTLxyUFFwMAiqukRNv+6f+9UAoAOFYoJw5LKHK6bI5W0iIi1zn0/Sdiu8u143QsCVk5+txoGSzfW8NN3leR4+T57sNIKREREZFCjJS6j9srpefOZIntnMwfAQBd+t+i6NzLr5vgljK5kvWb3d0rHxRpJ0ql6Nr45FgAgLmJTPrsrdFRW2EhcQCAKyIjRJopNB8AYMmTQ6V5ldKgp9Dab+lGg/zNPS5Q2j5j89/TJjTY6bL5wutH3sk6NU9TmpZHC0ZHvU/Wr5vlX9KGAgDeGDdfJJ0+uN7TRVJAY59SsFLaGEZKiYiIiBRipNR9WCklIiIiUspg0Dj6npPnN8ZtlVJr5+UlXz4n0jJLpBVz7sg+INIm3jIf/uD9aS+K7YfeeQQAcH1iJ72K02SondPT2nyXt0nuVnJtyzwAwLFo+SbTNTIQAHB5VO3KTjard0QFSP82ScHy4KaONt0BiDwtNTJB7yIQaeKoO19ujrx9qvgCAKBHnaP0YzAYYNBQwdRyTlPDSCkRERGRUlzRyW1cXim1Rq6s64l/kV0p9hVUSN8SukfJU0B886U0ZdTgG+e4uii6eWGyFB22vha+MFWSr3F2+qS2aWPF9qSiVQCAoir5vdozNhEA0GvYdABA9p9nxb6Uw5sAADMm3CPSfvthrVPlIXJG90H36l0EIpex/Zy86aZHncpry5dyK+agGx9s4EjlDEaDpnXstZzT1DBSSkRERKQQBzq5j9sqpaEB4QCAwc3lS5wul/qUNgsKEmlVlhp3FcHlrJP+B0cli7R2V/ZRlQfXNXcNZ1+/pJRWYvuulMb7NSe0TLLZvqfO/s79xjpVHiLybm98JE3tN3X8Ip1LQmq4Kjpqy2g0wahhInwt5zQ1jJQSERERKWUwaBtJz4FOjWKllIiIiEghNt+7j8srpdZm1ZquYwAAs2z+CBdKLwIAogPDRVpQYJhT1/PkIKK84jPSdRS+sazl8fc1zf+5egEAIMimE/eMCUv0Kg4Rkcv5S7O90s8jdjGrH6eEch9GSomIiIgUYqTUfdxWKW3WQlr3vdnAu0XaqWO/AwBqKstFWpvOXTXl7+jbnisGER3c+iEAYM2Jn0Xa78VlAIDBLaR106+J7agqT3//xplbWQEAKKux6FwSIiJyxPr5eP7IDyJtxo9fAgCKq6XfH+8gDwDt10+aDs/fP780YZ9St2G1nYiIiEghqfneqOGhrlK6detWjB49GsnJyTAYDFi7dq17npAXYaWUiIiISClNFVL1q0AVFxejW7duWLFihZueiPdR3Xw/+XV5zq+Z7a8AAPQYfJ+ic1MuU9fsrZYrmhmCw+MByE32APBDjtS2kRAsNX9M6jXM6ev4k6cmPQsA+H7jcp1LQkTkv+54bbbY3l8ozfG9f666++6ErV+K7RN7pa5X1SYpgvdg+Wmxb3u69Blouw49m/IlnhroNGLECIwYMUL1dXyZokqpxSL1FSwoKEBFqdwftKi4VKR7WkFB3T6lpgDn+zQWFhUDgN3zrCmT/vnLSw211/b88/UFxTYVeb5GnmF9na3/o+S9bO+jRFpU2nwuVZepez9ZPzOryypEmqV2PIClRvpsqymT7yOFhdLxpkD5eFd8xnojtffRwsJiaGlols6r+zcLDg5GcHCw6vz8kcGi4K9w+vRppKSkeKI8RKTBqVOn0KpVq8YPJN3wPkrk3Rq7j5aVlSE1NRVZWVmarxEREYGioiK7tEWLFmHx4sUNnmcwGLBmzRqMHTtW87V9gaJIaXJyMk6dOoXIyEjOs0XkRSwWCwoLC5GcnNz4waQr3keJvJPS+2hISAgyMzNRUVHR4HGNXevS/39GSWWKKqVGo5FRGCIvFR0drXcRSAHeR4m8l9L7aEhICEJCQtxcmqaLo++JiIiISHdc0YmIiIjIyxQVFeHo0aPi98zMTOzbtw9xcXFo3bq1jiVzH0UDnYiIiIjIc7799lsMHDiwTvqkSZOwatUqzxfIA1gpJSIiIiLdsU8pEREREemOlVIiIiIi0h0rpURERESkO1ZKiYiIiEh3rJQSERERke5YKSUiIiIi3bFSSkRERES6Y6WUiIiIiHTHSikRERER6Y6VUiIiIiLSHSulRERERKS7/w9H5TLTB5BxNgAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 800x300 with 3 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "{'axes': [<Axes: title={'center': '$T_s = 280\\\\,K$'}>,\n",
       "  <Axes: title={'center': '$T_s = 310\\\\,K$'}>],\n",
       " 'cbar_ax': <Axes: ylabel='mm/day'>}"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "make_masked_precip_panels(precip_data_dict=lpf_masked_prec,temp_cases=['280K', '310K'], time_steps=[46,-1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "dcfbc2ff-aa6c-4665-9564-144036d89e27",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Stored 'lpf_masked_prec' (dict)\n"
     ]
    }
   ],
   "source": [
    "%store lpf_masked_prec"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9df8fbc7-826c-46a6-947f-b195adee6650",
   "metadata": {},
   "source": [
    "## Fill holes function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "23707ac1-51cf-4554-9217-076c1a9ab583",
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.ndimage import binary_fill_holes\n",
    "from tqdm import tqdm\n",
    "from typing import Union, Dict\n",
    "\n",
    "def fill_holes(data: Union[Dict[str, Union[xr.DataArray, np.ndarray]], xr.DataArray, np.ndarray]) -> Union[Dict[str, xr.DataArray], xr.DataArray, np.ndarray]:\n",
    "    \"\"\"\n",
    "    Fill internal holes in 2D slices along the time axis.\n",
    "\n",
    "    Args:\n",
    "        data: A dict of DataArrays/arrays or a single DataArray/array (3D: time, y, x)\n",
    "\n",
    "    Returns:\n",
    "        Same structure as input with binary holes filled.\n",
    "    \"\"\"\n",
    "\n",
    "    def _fill_single_array(arr):\n",
    "        filled_slices = []\n",
    "        for t in range(arr.shape[0]):\n",
    "            filled = binary_fill_holes(arr[t]).astype(arr.dtype)\n",
    "            filled_slices.append(filled)\n",
    "        return np.stack(filled_slices, axis=0)\n",
    "\n",
    "    if isinstance(data, dict):\n",
    "        result = {}\n",
    "        for key, val in tqdm(data.items(), desc=\"Filling Holes\", unit=\"key\"):\n",
    "            arr = val.values if isinstance(val, xr.DataArray) else val\n",
    "            filled = _fill_single_array(arr)\n",
    "\n",
    "            if isinstance(val, xr.DataArray):\n",
    "                filled = xr.DataArray(filled, dims=val.dims, coords=val.coords, attrs=val.attrs)\n",
    "            result[key] = filled\n",
    "        return result\n",
    "\n",
    "    else:\n",
    "        arr = data.values if isinstance(data, xr.DataArray) else data\n",
    "        filled = _fill_single_array(arr)\n",
    "\n",
    "        if isinstance(data, xr.DataArray):\n",
    "            return xr.DataArray(filled, dims=data.dims, coords=data.coords, attrs=data.attrs)\n",
    "        return filled"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
